//这是一个自定义内核的图染色算法，主要针对的是HPCG中的线性方程组的系数矩阵A进行图染色。
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<cusparse.h>

//这个函数进行多哈希的处理
//我们多传入两个参数，一个是哈希函数的数量，一个数包含了所有取模的数组
//开的线程数量大于点的数量*哈希函数的数量
//一开始我们使用一个循环遍历所有的hash
//然后根据当前线程号选出适合自己的哈希，然后执行相同的计算就好了
__global__ void multi_hash_graph_coloring(int* A_device, size_t pitch,int* element_num_each_row_device, 
    int row_num, int* count_device, int color_num, int* coloring, int* old_coloring_device, int* hash_arr,
    int hash_arr_size){
    //计算当前节点编号
    int thread_id = threadIdx.x + blockIdx.x * blockDim.x;

    //每row_num个线程计算同一个哈希，所以首先根据当前的进程号做（thread_id/row_num）
    //得到需要使用的哈希函数取模的参数，然后根据参数进行计算
    
    int hash_num = thread_id / row_num;


    if(hash_num >= 0 && hash_num < hash_arr_size && thread_id < row_num * hash_arr_size){
        //这里说明是需要进行计算的
        int hash_parameters = hash_arr[hash_num];

        // if(hash_parameters == 23){
        //     printf("哈希值获取到23\n");
        // }

        //计算当前的点编号
        int point = thread_id % row_num;

        //如果周围的节点在上次没有染色，那么我们就进行染色
        if(old_coloring_device[point] == 0){
            //开始当前节点的染色工作

            //因为实际上开启的线程很多，所以我们需要计算一下当前线程对应的节点
            //我们做取模操作就可以得到当前节点所对应的点
            int col_num = element_num_each_row_device[point];

            //计算当前节点的哈希值，哈希函数是取模函数
            int random_num = (point + hash_parameters) % 29;

            //这个变量存储其他节点的哈希值
            int random_neighbor = 0;

            //用两个值判断当前函数在当前哈希下是极大值还是极小值
            int max = 1;
            int min = 1;

            //遍历所有的邻居节点
            int i;

            int neighbor_col;

            //和所有相邻节点相比较，极大值染色color + 1 + 2 * hash_num，极小值染色color + 2 + 2 * hash_num
            int* A_device_row = (int*)((char*)A_device + point * pitch);

            //遍历所有的邻居节点
            for(i = 0; i < col_num; i++){
                neighbor_col = A_device_row[i];

                //查看邻居节点有没有染色，如果已经染色了，就不要和邻居节点比
                if(old_coloring_device[neighbor_col] == 0){
                    //计算所有节点随机数
                    random_neighbor = (neighbor_col + hash_parameters) % 29;

                    //计算这个随机数和当前节点随机数的比较
                    if(random_neighbor > random_num && max == 1){
                        max = 0;
                    }

                    if(random_neighbor < random_num && min == 1){
                        min = 0;
                    }

                    //如果两个节点一样，如果都是极小值，那么节点编号比较小的就是比较小
                    //如果都是极大值，那么节点编号比较小的就是比较大
                    if(random_neighbor == random_num && point > neighbor_col){
                        min = 0;
                        max = 0;
                        break;
                    }
                }
            }


            //根据实际情况染色
            //如果一个节点又是极大又是极小说明他没有邻居了，按照极小处理
            if((min ==  1 || max == 1) && point < row_num && coloring[point] == 0){
                //看看需不需要染色
                // if(hash_parameters == 23){
                //     printf("point:%d,染色:%d\n", point);
                // }
                
                //我不知道这样子直接判断会不会有问题，但是这是唯一一个可以保证参数原子性更新的操作
                int old_color;

                if(max == 1){
                    old_color = atomicMax(&coloring[point], color_num + 1 + 2 * hash_num);
                    if(old_color == 0){
                        atomicAdd(count_device, 1);
                    }
                }else{
                    old_color = atomicMax(&coloring[point], color_num + 2 + 2 * hash_num);
                    if(old_color == 0){
                        atomicAdd(count_device, 1);
                    }
                }

                
            }
        }
    }
}

//检查颜色是不是打印正确的内核函数
//传入矩阵A以及每个节点的染色情况，看看某一个节点的邻居节点的颜色是不是冲突了
__global__ void check_graph_coloring(int* A_device, size_t pitch, int* element_num_each_row_device, 
    int row_num, int* coloring){
    //遍历每一行
    int thread_id = threadIdx.x + blockIdx.x * blockDim.x;
    
    if(thread_id < row_num){
        //将一行取出
        int* A_device_row = (int*)((char*)A_device + thread_id * pitch);

        int col_num = element_num_each_row_device[thread_id];

        int neighbor_col;

        int i;

        //获取当前节点的颜色
        int color = coloring[thread_id];

        for(i = 0; i < col_num; i++){
            //查看当前节点的列号
            neighbor_col = A_device_row[i];

            if(coloring[neighbor_col] != 0 && color != 0 && color == coloring[neighbor_col] && neighbor_col != thread_id){
                //这里说明染色冲突
                printf("出现错了，%d与%d冲突，冲突颜色%d\n", thread_id, neighbor_col, color);
            }
        }
    }
}

#define A_row_size 4096

int main(){
    //我们创造一个二维数组，这个二维数组一共4096行，27列，每一行就代表了图的节点数量
    //每一列就代表了就是一个点，一行中的不同元素代表了这个节点和哪些行有连接。
    //创造4096个指针
    int* A[A_row_size];

    //为每个指针创造27个空间
    //遍历变量
    printf("申请空间\n");

    int i;
    for(i = 0; i < A_row_size; i++){
        cudaMallocHost((void **)&A[i], 27 * sizeof(int));
    }

    //用一个数组来记录每一行元素的数量
    int element_num_each_row[A_row_size];

    //初始化每一行的数量
    for(i = 0 ; i < A_row_size; i ++){
        element_num_each_row[i] = 0;
    }
    
    //我们将数据写入。
    
    //从文件中读入数据
    FILE* Graph_File;

    //打开文件
    Graph_File = fopen("../data/A_Matrix_col_1.txt", "r");

    //用一个缓冲区来获取读取的内容
    char read_line_buf[1024];

    //用一个变量来遍历A的每一行
    int A_row_index = 0;

    printf("读取文件\n");

    fgets(read_line_buf, 1024, Graph_File);

    while (!feof(Graph_File)) {
        
        //用一个变量来定位A中每一行的每一个元素
        int A_col_index = 0;
        
        //分割字符串，这里是每一个子串
        char* part_char;

        part_char = strtok(read_line_buf, " ");


        while(part_char != NULL && A_col_index < 27){

            //第一个元素
            A[A_row_index][A_col_index] = atoi(part_char);

            // printf("%d ", A[A_row_index][A_col_index]);
            
            A_col_index++;
            element_num_each_row[A_row_index]++;

            part_char = strtok(NULL, " ");
        }

        //先读，才能查到eof
        fgets(read_line_buf, 1024, Graph_File);

        A_row_index++;
    }

    fclose(Graph_File);

    // printf("开始打印\n");
    // //打印A矩阵，看看情况
    
    // for(i = 0; i < 4096; i++){
    //     for(j = 0; j < element_num_each_row[i]; j++){
    //         printf("%d ", A[i][j]);
    //     }

    //     printf("\n");
    // }

    //这里打印一下element_num_each_row
    // for(i = 0; i < 4096; i++){
    //     printf("%d\n", element_num_each_row[i]);
    // }

    //开始进行图染色
    
    //首先在GPU端申请二维数组，使用cudaMallocPitch
    //https://www.cnblogs.com/csyisong/archive/2010/01/10/1643519.html
    
    //设备端A矩阵指针
    int* A_device = 0;
    //每一行对齐之后的大小
    size_t pitch;
    
    printf("申请设备端的A指针\n");
    //申请设备端的A指针
    cudaMallocPitch((void**)&A_device, &pitch, 27 * sizeof(int), A_row_size);

    //将数组A的数据拷贝给设备端，值的一提的是，数组A的每一行在存储上不是连续的，所以我们使用memcpy2D来处理
    //我们所能用的方式就是将A矩阵一行一行地传给设备端
    printf("pitch = %d\n", pitch);

    for(i = 0; i < A_row_size; i++){
        //因为pitch是字节数量，如果想让A_device这种int指针做地址偏移的话，需要将A_device的指针类型变化一下。
        cudaMemcpy((char *)A_device + i * pitch, A[i], 
            element_num_each_row[i] * sizeof(int), cudaMemcpyHostToDevice);
    }
    cudaThreadSynchronize();


    printf("将数组A的数据拷贝给设备端完毕\n");

    //然后将每一行多少个元素的数组拷贝到GPU端
    int* element_num_each_row_device = 0;

    cudaMalloc((void **)&element_num_each_row_device, A_row_size * sizeof(int));

    //将数据拷贝到设备端
    cudaMemcpy((void *)element_num_each_row_device, element_num_each_row, A_row_size * sizeof(int), cudaMemcpyHostToDevice);

    printf("然后将每一行多少个元素的数组拷贝到GPU端\n");

    //创建一个变量来保存现在染的数量
    int count = 0;

    int* count_device;

    //将这个变量传到显存上
    cudaMalloc((void**)&count_device, sizeof(int));

    cudaMemcpy(count_device, &count, sizeof(int), cudaMemcpyHostToDevice);

    //创建一个数组来表达每个节点所染的颜色
    int* coloring = 0;
    
    cudaMallocHost((void **)&coloring, sizeof(int) * A_row_size);

    int* coloring_device = 0;

    cudaMalloc((void **)&coloring_device, sizeof(int) * A_row_size);
    //将整块显存初始化为0
    cudaMemset(coloring_device, 0, sizeof(int) * A_row_size);

    //上次迭代的染色情况
    int* old_coloring_device = 0;

    cudaMalloc((void **)&old_coloring_device, sizeof(int) * A_row_size);

    //创建颜色数量
    int color_num = 0;

    //创建一个数组来存储哈希需要取模的值
    int* hash_arr_host = 0;
    int hash_arr_size = 2;

    cudaMallocHost((void **)&hash_arr_host, hash_arr_size * sizeof(int));

    //初始化两个哈希函数
    hash_arr_host[0] = 0;
    hash_arr_host[1] = 1;

    //将这个数组传到显存上
    int* hash_arr_device = 0;
    
    cudaMalloc((void **)&hash_arr_device, hash_arr_size * sizeof(int));

    cudaMemcpy(hash_arr_device, hash_arr_host, hash_arr_size * sizeof(int), cudaMemcpyHostToDevice);

    //创建线程结构，线性网格，每个block64个线程
    dim3 block_dimension(64);

    dim3 grid_dimension(A_row_size * hash_arr_size / 64 + 1);

    int old_count = -1;

    //循环调用内核函数
    while(count != A_row_size){
        printf("开始运行内核函数，count = %d\n", count);

        //将上次迭代的结果暂存起来
        cudaMemcpy(old_coloring_device, coloring_device, sizeof(int) * A_row_size, cudaMemcpyDeviceToDevice);

        old_count = count;

        multi_hash_graph_coloring<<<grid_dimension, block_dimension>>>(A_device, pitch, element_num_each_row_device, A_row_size, 
            count_device, color_num, coloring_device, old_coloring_device, hash_arr_device, hash_arr_size);

        cudaThreadSynchronize();

        color_num = color_num + 2 * hash_arr_size;
        
        //将count拷贝回来
        cudaMemcpy(&count, count_device, sizeof(int), cudaMemcpyDeviceToHost);
    }
    
    //查看错误
    //__global__ void check_graph_coloring(int* A_device, size_t pitch, int* element_num_each_row_device, 
    //int row_num, int* coloring)
    check_graph_coloring<<<grid_dimension, block_dimension>>>(A_device, pitch, element_num_each_row_device,
        A_row_size, coloring_device);

    cudaThreadSynchronize();
    
    printf("已染色颜色：%d，颜色种类：%d\n", count, color_num);

    //将颜色拷贝到cpu端
    cudaMemcpy(coloring, coloring_device, sizeof(int) * 4096, cudaMemcpyDeviceToHost);

    //打印所有的染色结果
    for(i = 0; i < 4096; i++){
        if(coloring[i] == 0){
            printf("%d:%d ", i, coloring[i]);
        }
    }

    // printf("\n");

    //析构A
    for(i = 0 ; i < A_row_size; i++){
        cudaFreeHost(A[i]);
    }

    //析构coloring_device
    cudaFree(coloring_device);

    //析构coloring
    cudaFreeHost(coloring);

    //析构element_num_each_row_device
    cudaFree(element_num_each_row_device);

    //析构设备端矩阵A
    cudaFree(A_device);

    //析构count_device
    cudaFree(count_device);

    //析构old_coloring_device
    cudaFree(old_coloring_device);

    //析构hash_arr_host
    cudaFreeHost(hash_arr_host);

    //析构hash_arr_device
    cudaFree(hash_arr_device);
}